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An efficient novel algorithm was developed to estimate the Density of States (DOS) for large systems by calculating the en- 
semble means of an extensive physical variable, such as the potential energy, U, in generalized canonical ensembles to inter- 
polate the interior reverse temperature curve /,(U) -20 , where S(U) is the logarithm of the DOS. This curve is computed 
with different accuracies in different energy regions to capture the dependence of the reverse temperature on U without setting 
prior grid in the U space. By combining with a U-compression transformation, we decrease the computational complexity from 
O(N*”) in the normal Wang Landau type method to O(N”) in the current algorithm, as the degrees of freedom of system N. 
The efficiency of the algorithm is demonstrated by applying to Lennard Jones fluids with various N, along with its ability to 
find different macroscopic states, including metastable states. 
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1 Introduction 


An essential problem for molecular simulations in complex 
systems with a rugged potential energy surface, such as 
proteins [1,2] and glasses [3] or first-order phase-coexist- 
ence systems, is how to adequately sample the conforma- 
tional space, since the existing high free-energy barriers 
prevent the regular Monte Carlo (MC) and Molecular Dy- 
namics (MD) simulations to visit all the important regions 
sufficiently. In attempting to overcome this difficulty, many 
enhanced sampling methods [4-18] have been developed; 
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most of them tried to modify the potential energy surface to 
make the distribution function in a preselected collective 
variable space be almost flat. Although these methods have 
achieved great accomplishments, more attempts are still 
needed to further improve the sampling efficiency, espe- 
cially in larger systems. 

The Density of States (DOS) of a system is defined as 


QU) = | dré(U —U(r)), where U(r) is the potential energy 


of the system (or more generally, any conformational func- 
tion). Here, r is a simple denotation of high-dimensional 
conformations, such as all-atomic position vectors, and the 
integral about r is over the whole conformational space. 
Usually, we define the microcanonical entropy S(U)=In 
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Q(U) and the interior reverse temperature £,(U)= a 


to describe the inner property of the system. Then, we have 
the flat-histogram ensemble, where the conformational dis- 
tribution Py(r xe “UO. The standard simulation techniques, 
such as MC or MD simulations, can generate the desired 
distribution if a new potential Vjew(r)=Ty S(U) is applied, 
where Ty is the temperature of thermostat, which could be 
any value. Thus the histogram in the U space P(U)xe*” 


Í drP(r)d(U—U(r)) is flat. Since the flat-histogram en- 


semble can visit a much wider energy range than the usual 
Canonical Ensemble (CE) where the conformational distri- 
bution P.(ryce PO , it has been widely used to overcome 
free-energy barriers and thus enhances the sampling for 
complex systems [4—6,9,10,13]. Here, @ =I/kgT, with kg 
Boltzmann constant and T the temperature. In the flat- 
histogram ensemble, S(U) is an input quantity and the target 
of simulations at the same time, so S(U) should be solved 
iteratively. The convergence efficiency and the accuracy of 
the iteration are essential for different flat-histogram meth- 
ods. 

Among them, the Wang Laudau (WL) algorithm [5] lifts 
the value of S in the current-visited U range at each simula- 
tion step until the visiting histogram in U space is almost 
flat. The WL algorithm is easy to implement and has been 
widely applied to various kinds of systems. However, the 
number of required simulation steps is at least in the order 
of O(N’) [8] with the system size N. The recently developed 
Statistical-Temperature Molecular Dynamics (STMD) 
method [10] updates the intensive variable T,(U) = 87" (U) 


instead of the extensive variable S(U) in the original WL, 
thus improving the efficiency of convergence to the order of 
O(N*”). A short point of these WL-like methods is that they 
do not satisfy the detailed balance condition thus with less 
numerical stability in applications. Alternatively, the multi- 
micro canonical method (MMCM) [7] calculates the mi- 
crocanonical temperature T„(E) and relates the T„(E) to the 
interior reverse temperature {,(E) based on the thermody- 
namic equivalence of differently defined temperatures and 
can converge within O(N*”) simulation steps. Different 
from the original WL algorithm and the STMD, the MMCM 
satisfies the detailed balance condition, then has some bene- 
fits in numerical stability. In the WL-like algorithms and the 
MMCM, it is necessary to cut the extensive variable U into 
homogeneous bins with a presumed bin size. We collect 
sufficient simulation data in each bin of U to estimate the 
value of £,(U) at the centers of each U bin, that is, preset a 
uniform grid {U;}, i=1,..., M and calculate £,(U) on the grid, 
then interpolate the whole curve £,(U) (or T,(U)). Since the 
varying rates of (U) at different U regions may be very 
different, the presetting of a suitable U grid for achieving 
both the accuracy and the efficiency is not trivial. We 
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should use a sparse {U;} grid in the slow-varied regions of 
£,(U), but a more dense grid in the fast-varied regions of 
(U). 

In this paper, we present the Fast Adaptive Flat-histo- 
gram Ensemble (FAFE) method to efficiently estimate 
(U). The FAFE calculates the ensemble means and fluctu- 
ations of U in a series of generalized ensembles to adap- 
tively form a U grid and to get the values of (U) on the 
grid, {(G, U}. The {U;} grid is non-uniform and the grid 
bin widths depend on both U and N automatically: dense 
{U;} points are applied in the £,(U) fast-changed regions but 
sparse {U;} in the U regions where Ø, varies slowly as U. 
In the U regions where £,(U) was already well estimated, 
we apply a compression transformation to further depress 
the visit time in that regions so that more simulation costs 
are focus on the others. By doing so, in the FAFE, the esti- 
mate of £,(U) can be converged within the order of N? 
simulation steps. Also, the FAFE satisfies the detailed bal- 
ance condition, and it can visit and identify coexistent 
phases, including metastable phases of systems. 


2 Theory and method 


When an estimate of the entropy S(U) is obtained, the CE 


means of any conformational function Q(r) can be easily 
calculated by 


O.w 

(O), = Lalm (1) 
> k Wi 

where Q;=0(r;) is the value of the Q(r) at the sampled con- 

formation rų under the flat-histogram simulation, and the 


ŠU (r, ))-U (m, kyl 


weight factor is w, =e . It is worth mention- 


ing that the weighted average is efficient only while the 
effective size of the weighted sample is large enough, which 
usually requires the generated sample to cover the target 
sample, that is, the visited U range of the flat-histogram 
ensemble sample covers that of the CE. We can approxi- 
mately estimate the statistical error of the weighted average 
by the effective size of the weighted sample: 


(Xm) 
Dm 


which decreases with increasing fluctuation of normalized 
weights. While the weights are constant, the effective size 
of the sample, B, is the realistic size of the sample, as per 
our expectation. Only while all significant regions of the 
target distribution are adequately sampled, the generated 
sample can be weighted efficiently to reproduce the target 
one, that is, the effective size B of the weighted sample is 
much larger than unity. 


B= (2) 
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For any f(U) ensemble, that is, Pine", the histo- 
gram Hither) could be approximately written as: 


H,(U) « oxo] sar -E 6) 
z 


by Taylor expansion, if N is very large, and UccN is an ex- 
tensive variable. The maximal point of HU) is approxi- 
mately equal to the ensemble average of U under the f(U) 
ensemble, Ey ~ <U>,. Thus, we have 


BA(E,) = f'(E;), 


' -2 n" (4) 

BYE,)=—-0;7 + f"(E,), 
where the single prime and double prime mean the first de- 
rivative and the second derivative, respectively. of is the 
fluctuation of U in the f(U) ensemble. We calculate E; = 


<U>; and o; =(U-E,)) in a series of different KU) 


ensembles to generate points on the curve £,(U). The whole 
curve £,(U) can be interpolated from these points, and 


S(U) can be estimated by integrating (U). 


Different f(U) ensembles, such as the CE f(U) =2U and 
the microcanonical ensemble 


f,.U)= le -U)@(E-U)| (5) 


can be used to construct 2,(U). Here, / is number of degrees 
of freedom, the Heaviside step function @(x) =1 while x>0, 
and zero otherwise. In this paper, we use the Generalized 
Canonical Ensemble (GCE) [17]: 


KOE AU +ZaU-U,) © 


to generate points on £,(U). In thermodynamic stable re- 
gions, Ø'<0. Thus, the distributions in the CE (@=0) are 


Gaussian-like functions around E;. However, in thermody- 
namic unstable regions such as two phase-coexistence re- 
gions, which cannot be sufficiently visited in the CE, we 
can use the GCE with a> J! (E) to stabilize the phase- 


coexistence regions, then reconstruct the data points of £,(U) 
in the regions. 

The implementation of our algorithm is straightforward. 
Without losing generality, we choose U as the total potential 
energy of a physical system and define T,(U)=B, (U). The 
steps of our algorithm are as follows: 

(1) choose an initial Tọ (usually the highest temperature 
we desired), and then set 7,(U)=T > and estimate S(U) = 
UIT); 

(2) simulate a segment of trajectory in the current S(U) 
ensemble and generate an equilibrium sample; 


(3) choose a few f(U) functions and calculate U, and 


fluctuations o; from the generated sample in the step 2 
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based on eq. (1). 

(4) Interpolate T,(U) and linearly extrapolate it out of the 
current visited regions from the ending U; from the obtained 
points {U,,T,)} (and the derivative of TU), {T/= ied / 
o? }), then integrate S(U) from the T,(U); 

(5) Repeat the steps 2, 3, and 4 until the whole T,(U) and 
S(U) are generated. 

At the step 3, we may use CEs, for example, f(U)= UIT; 
with a series of {7;}. Here, 7; is chosen from the current 
T,(U) function: the corresponding U ¿> which satisfies the 


equation T=T,( U ;), must locate in the current visited re- 


gions and |U-—U;-\|~o;-1, where o; is the fluctuation of U at 
the temperature T;_;. Also, we apply the mean and fluctua- 
tion at T; when and only when the effective size of the 
weighted sample, B, is sufficiently large. Note that we can 
calculate the means and fluctuations at many T; from the 
same sample without new simulations. The computational 
cost for the weighted mean is almost negligible compared 
with the generation of samples. In this step, we can check 
on whether larger æ should be applied in applying the GCE 
to estimate U +- For example, we first choose a=0 to cal- 


culate (T;, U, ), then set A =T", U= U,, and use a posi- 
tive æ to calculate the means in the GCE ensemble, and 
finally compare the results to determine whether a larger œ 
should be used to depress the possible deviation between 
the extreme point and the average U. The points (4, U;) can 
be calculated using many samples generated in different 
segments of simulations in the step 3: 


pe 
> B 


where U! and B’ are the weighted mean of U and the cor- 
responding B factor from the jth sample, respectively. The 


(7) 


statistical error of U , which is measured by the total B = 
> Bi factor, can be suppressed by this way. 


Before all regions are sufficiently visited, T,(U) in some 
regions may already converge if the B factor of the corre- 


sponding U is large enough. For example, in normal 
physical systems, 7,(U) at high temperatures can be ob- 
tained before the low-temperature regions are sufficient 
visited. 

The necessary visiting time in the high-U regions, which 
is proportional to N even in the flat-histogram ensemble, 
can be further suppressed, so the simulation can focus its 
sampling more in the low-U region. We transform U to a 


new variable Y by a ge (2%: 

dU b 
pressed function g,(z) slowly approaches zero when z— œ% 
and is almost unity when z<0. Here, U, and b are parameters 
to be determined. One possible selection is 


l where the com- 
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8g) =— (8) 
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with the parameter c~2 or 3. The U distribution of Y-flat 
samples is P(U)~1/(U-U.) when U>U,. At the end of each 
segment of simulations, we reset U, to be the lower bound- 
ary of the clear region and choose b in such a way that the 
expected visiting probability above U. is comparable with 


that below U., and then use S(Y) to generate the Y-flat 


sample. Here, SY) = S(U)-1ng, K = } Although the 


Y-flat simulations can still visit the high-U region ergodi- 
cally, the majority of the simulation time is spent in the 
low-U region to refine T,(U). The compression transfor- 
mation significantly decreases the required computational 
cost of our algorithm in very large systems. 


3 Results 


To verify the applicability of our algorithm, we applied it to 
a truncated-and-shifted Lennard Jones (LJ) fluid (r, = 2.5 
oj) with the reduced density p = 0.8. U was selected as the 
total potential energy, and the DOS were calculated in the 
range between 7=0.5 and T,=100, corresponding to a very 
large range of U in large systems. The ensemble mean of U 
is regarded as exact if the error of U is smaller than 5% of 
the fluctuation of U (i.e., B > 400). The Langevin thermostat 
was adopted to keep the system temperature constant and 
conformations were sampled every 10 MD steps. 7,(U) was 
calculated and updated every 10° steps along with the pa- 
rameters for the compressive transformation. In our simula- 
tions, the distance between neighboring U; is much larger 
than in the WL-like methods [10] in high-T regions and 
large-N systems but automatically decreases in low-T re- 
gions and small N systems. The results are shown in Figure 


=] 
1. The T,(U/N), c, =(=) , and s(u) = S(u)/N quickly 
u 


converged in the large U regions, so the scaling relations 
were easily found. Here u=U/N. Figure 2 shows the effi- 
ciency of our approaches by the number of required MD 
steps tcN” with y=0.5. Here, all ensemble means are re- 
quired to reach the condition B > 400 so that their relative 
errors are smaller than 0.05. 

As shown in Figure 3, at the lower temperature region, 
the obtained T,(U) is found to have multiple branches cor- 
responding to the DOS of some separated parts in the con- 
formational space. For a segment of simulated trajectory, 
we may know the branch of 7,(U) the trajectory belongs to 
by calculating the corresponding (T, <U>r) points. In our 
simulations, the trajectories generated in the right branch of 
T,(U) (liquid phase) have the expected flat histograms in the 
region U>U,, but quickly transit into the low-U region (sol- 
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Figure 1 (Color online) The obtained T,, EP , and § are shown as 
u 


functions of U in the LJ systems with different size N. The number of 
requiring points {(U;, 4;)} increases as N increases. 
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Figure 2 (Color online) The required MD steps for getting the same 
accuracy of £,(U) with the size of system N. The accuracy is measured by 
the B factors in estimating ensemble averages of U (see the text). The red 
dash line is from tupa N”. 


id phase). Similarly, simulations in the left branch of T,(U) 
(solid phase) homogeneously visit the U<U2 region and 
transit into the liquid phase. The two branches have an 
overlap region, indicating that the potential energy of the 
liquid phase can be lower than that of the solid phase, and U 
definitely cannot serve as the order parameter for the liq- 
uid/solid phase transition. It is very hard to visit (and more 
difficult to identify) the low-energy liquid and high-energy 
solid conformations with the simulation techniques other 
than ours. Therefore, our approach, which is able to visit 
more parts of the configurational space and to identify them 
by forming their own DOS without a priorly defined order 
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Figure 3 (Color online) T,(U) in the supersaturated liquid and the solid 
phases, respectively, of the LJ system with N=108. The transition can 
easily happen when U<U, (from liquid to solid) and U>U, (from solid to 
liquid). 
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Figure 4 (Color online) The 7,(U) function in a LJ system with N=108 
and p=0.8 at the low-temperature region. Each branch of the function cor- 
responds to a particular macroscopic phase. The blue-filled squares cover 


=f 
the unstable liquid/solid coexistence region where c, -(%) , which 
u 
can be thought as the heat capacity in the iso-U ensemble, are negative, as 
shown in the bottom panel. The green down-triangles correspond to the 
complete FCC solid, which is a metastable state in the current system. 
Some phase transitions from metastable solids to the stable solid (the black 
circles) can be found at lower temperatures. When 7,—0, more complex 
phase behaviors are expected. 


parameter, provides a new approach for better understand- 
ing phase transition processes. 


4 Discussion 


It is possible to calculate the total S(U) and T,(U) from the 
obtained TU) branches, S(U)= In >) exp[S'(U)], the 
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corresponding T` (U) =>) a (DITUN, where 
SU) =Si4 i) T! (E)dE (9) 


and T,'(U) are the ith branch of the entropy and temperature, 
respectively. Here o,(U)= exp[S'(U)]/ >) ,explS'(U)]. In 


the calculation, the integration constants {S'o} must be de- 
termined from another simulation technique. However, it 
should be noted that even in the exact S(U) ensemble, free 
transition between states still cannot be guaranteed, because 
the application of S(U) does not change P(r) inside each 
iso-U subspace. This is a general limitation in all low- 
dimension flat-histogram ensemble methods. The algo- 
rithms with good stabilities (e.g., satisfying the detailed 
balance in a long time) are helpful for analyzing and over- 
coming the limitation. For example, our approach generates 
TU) and S(U) (with an arbitrary constant) in each visited 
subspace instead of the expected total T,(U) and S(U) for 
the whole conformational space. Therefore, although the 
equilibrium among subspaces is difficult to reach, it still 
provides many useful results and might construct the total 
S(U) by determining the relative constants of S(U) for each 
subspace. Also, we may use the generalized ensembles de- 
scribed by eq. (6) to find more metastable states. For exam- 
ple, besides the stable liquid and solid phases in the LJ sys- 
tem, we detected some metastable solid states with different 
crystalline structures as well as the unstable multiphase co- 
existence, which are shown in Figure 4. These complex 
solid phases appear because the system is fixed at the den- 
sity p=0.8, so the interparticle distance in the Face-Center- 
Cubic (FCC) structure does not match the equilibrium dis- 
tance of the LJ potential, ro=2". 

The flat-histogram ensembles based on multiple collec- 
tive variables are expected to better enhance the sampling 
than the ensembles based on a single variable [11,19]. The 
computation complexity of our algorithm is O(N), where d 
is the number of the collective variables and y =1/2, which 
is much faster than the WL-like with y=1.5 to 3. Recently, 
Vogel et al. [20] combined the original WL algorithm with 
the Replica Exchange Method (REM) to form a parallel WL 
sampling method, which makes application of the WL algo- 
rithm in large systems much efficient. It is also possible to 
combine the FAFE algorithm with the REM to further im- 
prove the efficiency of sampling. We may apply replicas in 
different U ranges [20], or in different branches of T,(U), 
then exchange conformations among these replicas. It might 
provide a possibility to overcome the main limitation in the 
application of REM, which requires O(N'”) replicas to 
achieve sufficient exchange acceptance ratio [21]. 
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